Figures

1 Import library

Code

# data viz
library(ggplot2)
library(ggh4x)
library(ggOceanMaps)
library(patchwork)
library(viridis)
library(ggpubr)
library(grid)
library(viridis)
library(ggnewscale)
library(gridExtra)
library(png)

# table
library(gt)
library(gtable)

# map
library(ggmap)
library(ggsn)
library(sf)
library(sp)
library(smoothr)

# kernel density
library(eks)
library(ggdensity)

# stat
library(Hmisc)
library(circular)
library(CircStats)

# char
library(stringr)

# solar angle and timezone
library(suncalc)
library(lutz)

# data wrangling
library(tidyr)
library(dplyr)
library(data.table)
library(purrr)
library(forcats)
library(lubridate)

2 Setting up custom function

2.1 windrose

Code

windrose <-
  function(data_to_plot,
           grid = NULL,
           set_title = NULL,
           legend_position = "none",
           facet_wrap = F,
           facet_grid = F,
           hour_sunset = 20,
           hour_sunrise = 7) {
    # this code comes from Roxanne
    uniqhours <- 1:24 * (360 / 24)

    # trick to align hours om the graph
    data_to_plot <- rbind(data_to_plot[2:nrow(data_to_plot), ], data_to_plot[1, ])

    for (i in 1:24) {
      # turn hours to radians
      if (i == 1) {
        temp <- rep(uniqhours[i], data_to_plot$nb_ind_hour[i])
        day_night <- rep("night", data_to_plot$nb_ind_hour[i])
      } else {
        temp <- c(temp, rep(uniqhours[i], data_to_plot$nb_ind_hour[i]))
        day_night <- c(day_night, rep(
          if_else(between(i, hour_sunrise, hour_sunset), "day", "night"),
          data_to_plot$nb_ind_hour[i]
        ))
      }
    }
    data2 <- data.frame(direction = temp)

    deg <- 15 # choose bin size (degrees/bin)
    dir.breaks <- seq(0 - (deg / 2), 360 + (deg / 2), deg) # define the range of each bin

    # assign each direction to a bin range
    dir.binned <-
      cut(data2$direction,
        breaks = dir.breaks,
        ordered_result = TRUE
      )
    # generate pretty lables
    dir.labels <- as.character(c(seq(0, 360 - deg, by = deg), 0))

    # replace ranges with pretty bin labels
    levels(dir.binned) <- dir.labels

    # Assign bin names to the original data set
    data2$dir.binned <- dir.binned

    # add origin if any
    if (facet_wrap == "origin") {
      data2$origin <- unique(data_to_plot$origin)
    } else if (facet_wrap == "trip") {
      data2$trip <- unique(data_to_plot$trip)
    }

    # add trip and origin if any
    if (facet_grid) {
      data2$origin <- unique(data_to_plot$origin)
      data2$trip <- unique(data_to_plot$trip)
    }

    # set up max value
    maxvalue <- 20

    # initialise the plot
    plt.dirrose_2 <- ggplot()

    # check if grid
    if (!is.null(grid)) {
      plt.dirrose_2 <- plt.dirrose_2 +
        geom_hline(
          yintercept = grid,
          colour = "grey20",
          linewidth = .2
        )
    }
    plt.dirrose_2 <- plt.dirrose_2 +
      geom_vline(
        xintercept = c(seq(1, 24, 2)),
        colour = "grey30",
        linewidth = 0.2
      ) + # 24 vertical lines at center of the 30? ranges.
      geom_hline(
        yintercept = maxvalue,
        colour = "black",
        linewidth = .5
      ) + # Darker horizontal line as the top border (max).
      # On top of everything we place the histogram bars.
      geom_bar(
        data = data2,
        aes(x = dir.binned, fill = day_night),
        width = 1,
        colour = "black",
        linewidth = 0.3
      ) +
      # geom_bar(data = data2, aes(x = dir.binned2), width = 1, colour="black", size = 0.3,fill="salmon",alpha=0.9) +
      scale_x_discrete(
        drop = FALSE,
        labels = c(
          0, "", 2, "", 4, "", 6, "", 8, "", 10, "", 12, "", 14, "", 16, "", 18, "", 20, "", 22, ""
        )
      ) +
      scale_fill_manual(values = c("white", "darkgrey"), name = "Time of day") +
      labs(x = "Time (hours)", y = "Count", title = set_title) +
      coord_polar(start = -(deg / 2) * (pi / 180))

    # if facet
    if (facet_wrap == "origin") {
      plt.dirrose_2 <- plt.dirrose_2 +
        facet_wrap2(. ~ origin)
    } else if (facet_wrap == "trip") {
      plt.dirrose_2 <- plt.dirrose_2 +
        facet_wrap2(trip ~ .,
          strip.position = "right"
        )
    }

    # if facet
    if (facet_grid) {
      plt.dirrose_2 <- plt.dirrose_2 +
        facet_grid2(trip ~ origin)
    }

    # Wraps the histogram into a windrose
    plt.dirrose_2 <- plt.dirrose_2 +
      theme_bw() +
      theme(
        legend.position = legend_position,
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        panel.background = element_blank(),
        legend.key.size = unit(0.5, "cm")
      )

    # return
    return(plt.dirrose_2)
  }

2.2 matlab_to_posix

Code

# function to convert matlab date into R posix
matlab_to_posix <- function(x, timez = "UTC") {
  days <- x - 719529 # 719529 = days from 1-1-0000 to 1-1-1970
  secs <- days * 86400 # 86400 seconds in a day

  return(as.POSIXct(
    strftime(
      as.POSIXct(secs,
        origin = "1970-1-1",
        tz = "UTC"
      ),
      format = "%Y-%m-%d %H:%M:%S",
      tz = "UTC",
      usetz = FALSE
    ),
    tz = timez
  ))
}

3 Import data

Let’s load data_dive, i.e. the output of data_wrangling.qmd, and filter only on animals leaving from Ano Nuevo.

Code

# import the data
data_dive <- readRDS("../export/data_dive.rds")

# filter on seals departing from ANNU
data_dive <- data_dive %>%
  filter(DepartureLocation == "ANNU")

# remove id == `2016014` and id = `2018002` (no benthic dives)
data_dive <- data_dive %>%
  filter(!id %in% c("2016014", "2018002"))

Code

# change few things
data_dive <- data_dive %>%
  mutate(
    # change divetype
    DiveType = data.table::fifelse(DiveType == 3 &
      (abs(bathy) - Maxdepth) > 100, 0, DiveType),
    DiveTypeName = data.table::fifelse(
      DiveTypeName == "Benthic" &
        (abs(bathy) - Maxdepth) > 100,
      "Transit",
      DiveTypeName
    )
  )

# add Post Molting vs Post Breeding trip
data_dive <- data_dive %>%
  group_by(id) %>%
  mutate(trip = if_else(between(month(first(date)), 1, 3), "Post Breeding", "Post Molting"))

4 Data Visualisation

4.1 Table 1

We need to define a custom function to calculate the summarize version of each column differently, i.e. weighted mean and sum.

Code

compute_table_1_a <- function(x) {
  # if the sum of the column is below 100
  if (sum(x) < 100) {
    # it's definitely a proportion
    # so we weighted average these
    y <- weighted.mean(x, c(1068746, 2207361))
    y <- paste0(round(y * 100, 1), "%")
    return(y)
    # otherwise we just sum the column
  } else {
    y <- prettyNum(sum(x), big.mark = ",")
    return(y)
  }
}

Code

# shark and orca area = NA if no location data
data_dive <-
  data_dive %>%
  mutate(
    shark_area = if_else(is.na(Lat), NA, shark_area),
    orca_area = if_else(is.na(Lat), NA, orca_area)
  ) %>%
  # add year of deployment
  group_by(id) %>%
  mutate(
    year = first(format(date, format = "%Y")),
    shark_and_orca = if_else(is.na(Lat),
      NA,
      if_else(shark_area + orca_area == 2, 1, 0)
    ),
    shark_or_orca = if_else(is.na(Lat),
      NA,
      if_else(shark_area + orca_area >= 1, 1, 0)
    )
  )

# build table
table_1_a <- data_dive %>%
  group_by(trip, id) %>%
  summarise(
    n_dives_ind_all = length(DiveNumber),
    n_dives_ind_w_loc = length(DiveNumber[!is.na(Lat)]),
    n_dives_ind_wo_loc = length(DiveNumber[is.na(Lat)]),
    n_benthic_dives_all = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic"])
    ),
    n_benthic_dives_w_loc = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_in_shark = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        shark_area == 1 &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_in_orca = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        orca_area == 1 &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_shark_and_orca = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        shark_and_orca == 1 &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_shark_or_orca = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        shark_or_orca == 1 &
        !is.na(Lat)])
    ),
    .groups = "drop"
  ) %>%
  group_by(trip) %>%
  summarise(
    n_ind = length(id),
    n_dives_all = sum(n_dives_ind_all, na.rm = T),
    n_dives_w_loc = sum(n_dives_ind_w_loc, na.rm = T),
    n_dives_wo_loc = sum(n_dives_ind_wo_loc, na.rm = T),
    n_benthic_dives = sum(n_benthic_dives_all, na.rm = T),
    n_benthic_dives_proportion = sum(n_benthic_dives, na.rm = T) / sum(n_dives_all, na.rm = T),
    n_benthic_dives_in_shark_proportion = sum(n_benthic_dives_ind_in_shark, na.rm = T) /
      sum(n_benthic_dives_w_loc, na.rm = T),
    n_benthic_dives_in_orca_proportion = sum(n_benthic_dives_ind_in_orca, na.rm = T) /
      sum(n_benthic_dives_w_loc, na.rm = T),
    n_benthic_dives_in_shark_and_orca = sum(n_benthic_dives_ind_shark_and_orca, na.rm = T) /
      sum(n_benthic_dives_w_loc, na.rm = T),
    n_benthic_dives_in_shark_or_orca = sum(n_benthic_dives_ind_shark_or_orca, na.rm = T) /
      sum(n_benthic_dives_w_loc, na.rm = T)
  ) %>%
  gt(rowname_col = "trip") %>%
  tab_spanner(
    label = "All dives",
    columns = c(
      "n_dives_all",
      "n_dives_w_loc",
      "n_dives_wo_loc"
    )
  ) %>%
  tab_spanner(
    label = "Benthic dives",
    columns = c(
      "n_benthic_dives",
      "n_benthic_dives_proportion"
    )
  ) %>%
  tab_spanner(
    label = "% Benthic dives per area",
    columns = c(
      "n_benthic_dives_in_shark_proportion",
      "n_benthic_dives_in_orca_proportion",
      "n_benthic_dives_in_shark_or_orca",
      "n_benthic_dives_in_shark_and_orca"
    )
  ) %>%
  # rename columns
  cols_label(
    trip = "Trip",
    n_ind = "# seals",
    n_dives_all = "Total",
    n_dives_w_loc = "w/ loc",
    n_dives_wo_loc = "w/o loc",
    n_benthic_dives = "Total",
    n_benthic_dives_proportion = "Proportion",
    n_benthic_dives_in_shark_proportion = "Shark",
    n_benthic_dives_in_orca_proportion = "Orca",
    n_benthic_dives_in_shark_or_orca = "Total",
    n_benthic_dives_in_shark_and_orca = "Overlap"
  ) %>%
  # summary rows
  grand_summary_rows(
    columns = everything(),
    fns = list(label = "OVERALL", fn = "compute_table_1_a")
  ) %>%
  fmt_percent(
    columns = c(
      "n_benthic_dives_in_shark_proportion",
      "n_benthic_dives_in_orca_proportion",
      "n_benthic_dives_proportion",
      "n_benthic_dives_in_shark_or_orca",
      "n_benthic_dives_in_shark_and_orca"
    ),
    decimal = 1
  ) %>%
  fmt_number(
    columns = c(
      "n_dives_all",
      "n_dives_w_loc",
      "n_dives_wo_loc",
      "n_benthic_dives"
    ),
    decimal = 0
  ) %>%
  # color rows
  opt_row_striping() %>%
  # width
  tab_options(table.width = pct(130))

Code

# custom function to summarize each column differently
compute_table_1_b <- function(x) {
  # if the column is numeric
  if (all(is.numeric(x))) {
    # and the sum is below 100
    if (sum(x) < 100) {
      # it is time
      # so we "circular" weighted average them
      y <- psych::circadian.mean(
        c(
          rep(x[1], 237),
          rep(x[2], 187)
        ),
        na.rm = T
      )
      y <- round(y, 1)
    } else {
      # otherwise it is just number
      # so we "just" weighted average them
      y <- weighted.mean(x, c(237, 187))
      y <- round(y, 1)
    }
  } else {
    y <- "\u00b1"
  }
  return(y)
}

Code

# shark and orca area = NA if no location data
data_dive <- data_dive %>%
  mutate(
    shark_area = if_else(is.na(Lat), NA, shark_area),
    orca_area = if_else(is.na(Lat), NA, orca_area)
  ) %>%
  # add year of deployment
  group_by(id) %>%
  mutate(
    year = first(format(date, format = "%Y")),
    shark_and_orca = if_else(is.na(Lat),
      NA,
      if_else(shark_area + orca_area == 2, 1, 0)
    ),
    shark_or_orca = if_else(is.na(Lat),
      NA,
      if_else(shark_area + orca_area >= 1, 1, 0)
    )
  )

# build table
table_1_b <- data_dive %>%
  group_by(trip, id) %>%
  summarise(
    n_dives_ind_all = length(DiveNumber),
    n_dives_ind_w_loc = length(DiveNumber[!is.na(Lat)]),
    n_dives_ind_wo_loc = length(DiveNumber[is.na(Lat)]),
    n_benthic_dives_all = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic"])
    ),
    n_benthic_dives_w_loc = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_in_shark = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        shark_area == 1 &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_in_orca = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        orca_area == 1 &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_shark_and_orca = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        shark_and_orca == 1 &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_shark_or_orca = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        shark_or_orca == 1 &
        !is.na(Lat)])
    ),
    hour_day_departure = first(hour(date_tz) +
      minute(date_tz) / 60 +
      second(date_tz) / (60 * 60)),
    hour_day_arrival = last(hour(date_tz) +
      minute(date_tz) / 60 +
      second(date_tz) / (60 * 60)),
    .groups = "drop"
  ) %>%
  group_by(trip) %>%
  summarise(
    n_dives_mean = mean(n_dives_ind_all, na.rm = T),
    n_dives_plus_minus = "\u00b1",
    n_dives_sd = sd(n_dives_ind_all, na.rm = T),
    n_benthic_dives_mean = mean(n_benthic_dives_all, na.rm = T),
    n_benthic_dives_plus_minus = "\u00b1",
    n_benthic_dives_sd = sd(n_benthic_dives_all, na.rm = T),
    hour_day_departure_mean = psych::circadian.mean(hour_day_departure, na.rm = T),
    hour_day_departure_plus_minus = "\u00b1",
    hour_day_departure_sd = psych::circadian.sd(hour_day_departure, na.rm = T)$sd,
    hour_day_arrival_mean = psych::circadian.mean(hour_day_arrival, na.rm = T),
    hour_day_arrival_plus_minus = "\u00b1",
    hour_day_arrival_sd = psych::circadian.sd(hour_day_arrival, na.rm = T)$sd
  ) %>%
  gt(rowname_col = "trip") %>%
  tab_spanner(
    label = "# dives",
    columns = c(
      "n_dives_mean",
      "n_dives_plus_minus",
      "n_dives_sd"
    )
  ) %>%
  tab_spanner(
    label = "# benthic dives",
    columns = c(
      "n_benthic_dives_mean",
      "n_benthic_dives_plus_minus",
      "n_benthic_dives_sd"
    )
  ) %>%
  tab_spanner(
    label = "Hour departure",
    columns = c(
      "hour_day_departure_mean",
      "hour_day_departure_plus_minus",
      "hour_day_departure_sd"
    )
  ) %>%
  tab_spanner(
    label = "Hour arrival",
    columns = c(
      "hour_day_arrival_mean",
      "hour_day_arrival_plus_minus",
      "hour_day_arrival_sd"
    )
  ) %>%
  # tab_spanner(
  #   label = "per individual",
  #   columns = c(
  #     "n_dives_mean",
  #     "n_dives_plus_minus",
  #     "n_dives_sd",
  #     "n_benthic_dives_mean",
  #     "n_benthic_dives_plus_minus",
  #     "n_benthic_dives_sd",
  #     "hour_day_departure_mean",
  #     "hour_day_departure_plus_minus",
  #     "hour_day_departure_sd",
  #     "hour_day_arrival_mean",
  #     "hour_day_arrival_plus_minus",
  #     "hour_day_arrival_sd"
  #   )
  # )  %>%
  # rename columns
  cols_label(
    trip = "Trip",
    n_dives_mean = md("Mean"),
    n_dives_plus_minus = "\u00b1",
    n_dives_sd = md("SD"),
    n_benthic_dives_mean = md("Mean"),
    n_benthic_dives_plus_minus = "\u00b1",
    n_benthic_dives_sd = md("SD"),
    hour_day_departure_mean = md("Mean"),
    hour_day_departure_plus_minus = "\u00b1",
    hour_day_departure_sd = md("SD"),
    hour_day_arrival_mean = md("Mean"),
    hour_day_arrival_plus_minus = "\u00b1",
    hour_day_arrival_sd = md("SD")
  ) %>%
  fmt_number(
    columns = c(
      "n_dives_mean",
      "n_dives_sd",
      "n_benthic_dives_mean",
      "n_benthic_dives_sd",
      "hour_day_departure_mean",
      "hour_day_departure_sd",
      "hour_day_arrival_mean",
      "hour_day_arrival_sd"
    ),
    decimal = 1,
    use_seps = FALSE
  ) %>%
  # alignement
  cols_align(
    columns = c(
      "n_dives_sd",
      "n_benthic_dives_sd",
      "hour_day_departure_sd",
      "hour_day_arrival_sd"
    ),
    align = "left"
  ) %>%
  # summary rows
  grand_summary_rows(
    columns = everything(),
    fns = list(label = "OVERALL", fn = "compute_table_1_b")
  ) %>%
  # color rows
  opt_row_striping()

Code

# display first table
table_1_a

Table 1:

Summary of dataset characteristics for Post Breeding and Post Molting foraging trips in northern elephant seals.

	# seals	All dives			Benthic dives		% Benthic dives per area
		Total	w/ loc	w/o loc	Total	Proportion	Shark	Orca	Total	Overlap
Post Breeding	235	1,058,860	1,032,565	26,295	36,443	3.4%	79.8%	100.0%	100.0%	79.8%
Post Molting	187	2,207,361	2,032,528	174,833	44,594	2.0%	90.5%	100.0%	100.0%	90.5%
OVERALL	422	3,266,221	3,065,093	201,128	81,037	2.5%	87%	100%	100%	87%

Code

# display second table
table_1_b

Table 2:

Individual-level summary of dataset characteristics for Post Breeding and Post Molting foraging trips in northern elephant seals

	# dives			# benthic dives			Hour departure			Hour arrival
	Mean	±	SD	Mean	±	SD	Mean	±	SD	Mean	±	SD
Post Breeding	4505.8	±	905.7	155.1	±	275.7	19.0	±	1.8	15.8	±	2.5
Post Molting	11804.1	±	1379.4	238.5	±	471.3	19.2	±	1.7	13.8	±	2.0
OVERALL	7724.6	±	1114.6	191.9	±	362	19.1	±	1.8	14.9	±	2.3

Notes

• Table 1: Summary of dataset characteristics for Post Breeding and Post Molting foraging trips in northern elephant seals.

• Table 2: Individual-level summary of dataset characteristics for Post Breeding and Post Molting foraging trips in northern elephant seals

Code

# get the hours of departure and arrival
hour_data <- data_dive %>%
  group_by(trip, id) %>%
  summarise(
    hour_day_departure = first(hour(date_tz) +
      minute(date_tz) / 60 +
      second(date_tz) / (60 * 60)),
    hour_day_arrival = last(hour(date_tz) +
      minute(date_tz) / 60 +
      second(date_tz) / (60 * 60)),
    .groups = "drop"
  ) %>%
  filter(!is.na(hour_day_departure) & !is.na(hour_day_arrival))

# compare the circular mean hour of departure between two types of trip
watson.williams.test(
  list(
    pb = hour_data %>%
      filter(trip == "Post Breeding") %>%
      pull(hour_day_departure) %>%
      circular(., units = "hours"),
    pm = hour_data %>%
      filter(trip == "Post Molting") %>%
      pull(hour_day_departure) %>%
      circular(., units = "hours")
  )
)

    Watson-Williams test for homogeneity of means

data:  pb and pm
F = 0.12556, df1 = 1, df2 = 399, p-value = 0.7233
sample estimates:
Circular Data: 
Type = angles 
Units = hours 
Template = none 
Modulo = asis 
Zero = 0 
Rotation = counter 
mean of pb mean of pm 
 -5.033036  -4.759905 

Code

watson.wheeler.test(
  list(
    pb = hour_data %>%
      filter(trip == "Post Breeding") %>%
      pull(hour_day_departure) %>%
      circular(., units = "hours"),
    pm = hour_data %>%
      filter(trip == "Post Molting") %>%
      pull(hour_day_departure) %>%
      circular(., units = "hours")
  )
)

    Watson-Wheeler test for homogeneity of angles

data:  pb and pm
W = 0.14773, df = 2, p-value = 0.9288

More here: https://stats.stackexchange.com/questions/265865/comparisons-of-circular-means

4.2 Figure 1

Individual are not chosen completely randomly, we picked randomly individuals that performed ~180 dives per 36h.

Code

# to choose an individual ~180 dives for 36 hours
data_dive %>%
  # rename date as datetime
  mutate(datetime = date) %>%
  # then by individual
  group_by(id) %>%
  # calculate
  mutate( # the difference time with the first dive
    diff_start = abs(as.numeric(difftime(
      first(datetime), datetime,
      units = "days"
    ))),
    # the difference time with the last dive
    diff_end = abs(as.numeric(difftime(
      last(datetime), datetime,
      units = "days"
    )))
  ) %>%
  filter(diff_start <= nb_days | diff_end <= nb_days) %>%
  mutate(phase = factor(
    if_else(diff_start <= nb_days, "Departure", "Arrival"),
    levels = c("Departure", "Arrival")
  )) %>%
  group_by(id, trip) %>%
  summarise(nb_dives = length(DiveNumber)) %>%
  arrange(nb_dives) %>%
  View()

Code

# import the csv export
tdr_data <- lapply(
  list.files(
    path = "../data/",
    pattern = "*_iknos_raw_data.csv",
    full.names = T
  ),
  function(x) {
    # extract id
    id_name <- last(strsplit(strsplit(x, split = "_")[[1]][1],
      split = "/"
    )[[1]])
    # load file
    dt <- readr::read_csv(x)
    # add id
    dt %>%
      mutate(id = as.numeric(id_name))
  }
) %>% bind_rows()

# choose the lag
nb_days <- 1.5

# compute the data
data_plot <- tdr_data %>%
  # convert time into date, a Posix
  mutate(
    datetime = matlab_to_posix(time),
    date = as.Date(datetime)
  ) %>%
  # add beginning and end date of the trip
  merge(
    .,
    data_dive %>%
      group_by(id, trip) %>%
      summarise(
        date_low = min(date),
        date_high = max(date),
        .groups = "drop"
      ),
    by = "id",
    all.x = TRUE
  ) %>%
  # identify what to keep
  mutate(to_keep = if_else(between(datetime, date_low, date_high), T, F)) %>%
  # keep only between the beginning and the end
  filter(to_keep) %>%
  # sort to make sure the difftime are being done correctly
  arrange(id, time) %>%
  # merge to get lat and lon
  merge(
    .,
    data_dive %>%
      group_by(id, date = as.Date(date)) %>%
      summarise(
        lat = median(Lat, na.rm = T),
        lon = median(Long, na.rm = T),
        .groups = "drop"
      ),
    by = c("id", "date"),
    all.x = TRUE
  ) %>%
  # then by individual
  group_by(id) %>%
  # calculate
  mutate( # the difference time with the first dive
    diff_start = abs(as.numeric(difftime(
      first(datetime), datetime,
      units = "days"
    ))),
    # the difference time with the last dive
    diff_end = abs(as.numeric(difftime(
      last(datetime), datetime,
      units = "days"
    )))
  ) %>%
  filter(diff_start <= nb_days | diff_end <= nb_days) %>%
  mutate(phase = factor(
    if_else(diff_start <= nb_days, "Departure", "Arrival"),
    levels = c("Departure", "Arrival")
  ))

# add bathy
setDT(data_plot)
setDT(data_dive)
data_dive[, datetime := date]
setkeyv(data_plot, c("id", "datetime"))
setkeyv(data_dive, c("id", "datetime"))
data_plot <- copy(data_dive) %>%
  .[, .(id, datetime, bathy, DiveTypeName)] %>%
  .[data_plot, roll = T]

# get sunrise sunset information
sun_data <- data_plot %>%
  select(date, lat, lon, id, phase) %>%
  unique() %>%
  ungroup() %>%
  mutate(getSunlightTimes(data = ., keep = c("sunrise", "sunset"))) %>%
  # fuck it, getSunlightTimes gives the sunset of next day... so we subtract 1
  mutate(sunset = if_else( # if sunset is the next day
    abs(as.numeric(
      difftime(as.POSIXct(date), sunset, units = "days")
    )) > 1,
    # subtract 1
    sunset - (1 * 60 * 60 * 24),
    # otherwise leave it that way
    sunset
  ))

# trim if night beging before the trip or end after they arrive
sun_data <-
  merge(
    sun_data %>% select(date, lat, lon, id, phase, sunrise, sunset),
    data_plot %>%
      group_by(id, phase, trip) %>%
      summarise(
        date_low = min(datetime),
        date_high = max(datetime),
        .groups = "drop"
      ),
    by = c("id", "phase"),
    all.x = TRUE,
  ) %>%
  mutate(
    # make sure sunset is not outside the 36h period
    correct_sunset = if_else(
      # if not
      between(sunset, date_low, date_high),
      sunset,
      # if so, replace by one extremum
      if_else(sunset < date_low, date_low, date_high)
    ),
    # make sure sunrise is not outside the 36h period
    correct_sunrise = if_else(
      # if not
      between(sunrise, date_low, date_high),
      sunrise,
      # if so, replace by one extremum
      if_else(sunrise < date_low, date_low, date_high)
    )
  ) %>%
  # shift everything to start 2000-01-01 00:00:00
  mutate(
    correct_sunset_shift = correct_sunset - (
      date_low - as.POSIXct("2000-01-01 00:00:00", tz = "UTC")
    ),
    correct_sunrise_shift = correct_sunrise - (
      date_low - as.POSIXct("2000-01-01 00:00:00", tz = "UTC")
    )
  )

# plot
departure_dp <- data_plot %>%
  filter(phase == "Departure") %>%
  group_by(id, phase) %>%
  # make the date starts on 2000-01-01 00:00:00
  mutate(datetime = datetime - (
    first(datetime) - as.POSIXct("2000-01-01 00:00:00", tz = "UTC")
  )) %>%
  ungroup() %>%
  mutate(`Dive Type` = if_else(DiveTypeName == "Benthic", "Benthic", "Other")) %>%
  ggplot(aes(x = datetime, y = depth, col = `Dive Type`, group = id)) +
  scale_color_manual(values = c("Benthic" = "#008c49", "Other" = "grey10")) +
  geom_rect(
    data = sun_data %>%
      filter(phase == "Departure"),
    aes(
      xmin = correct_sunset_shift,
      xmax = correct_sunrise_shift,
      ymin = -Inf,
      ymax = +Inf
    ),
    alpha = 0.4,
    inherit.aes = FALSE
  ) +
  # geom_path(aes(x = datetime, y = abs(bathy))) +
  geom_path() +
  scale_x_datetime(
    # date_breaks = "6 hours",
    # date_labels = "+ %Hh",
    breaks = seq.POSIXt(
      as.POSIXct("2000-01-01 06:00:00", tz = "UTC"),
      as.POSIXct("2000-01-02 06:00:00", tz = "UTC"),
      by = "6 hour"
    ),
    labels = c("+6h", "+12h", "+18h", "+24h", "+30h"),
    position = "top"
  ) +
  scale_y_reverse(limits = rev(range(data_plot$depth))) +
  labs(x = "Time", y = "Depth (m)") +
  # facet_nested(trip + id ~ phase,
  facet_nested(id ~ phase,
    scales = "free_x",
    # independent = "x"
  ) +
  theme(
    strip.placement = "outside",
    strip.text.y = element_blank(),
    legend.position = "top",
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank(),
    panel.grid.major = element_line(color = "grey90"),
    axis.line.x = element_line(
      color = "black",
      # arrow = arrow(length = unit(0.2, "lines"), type = "closed")
    ),
    axis.line.y = element_line(
      color = "black",
      arrow = arrow(
        length = unit(0.2, "lines"),
        type = "closed",
        ends = "first"
      )
    ),
    axis.title.x = element_blank()
  )

arrival_dp <- data_plot %>%
  filter(phase == "Arrival") %>%
  group_by(id, phase) %>%
  # make the date starts on 2000-01-01 00:00:00
  mutate(datetime = datetime - (
    first(datetime) - as.POSIXct("2000-01-01 00:00:00", tz = "UTC")
  )) %>%
  ungroup() %>%
  mutate(`Dive Type` = if_else(DiveTypeName == "Benthic", "Benthic", "Other")) %>%
  ggplot(aes(x = datetime, y = depth, col = `Dive Type`, group = id)) +
  scale_color_manual(values = c("Benthic" = "#008c49", "Other" = "grey10")) +
  geom_rect(
    data = sun_data %>%
      filter(phase == "Arrival"),
    aes(
      xmin = correct_sunset_shift,
      xmax = correct_sunrise_shift,
      ymin = -Inf,
      ymax = +Inf
    ),
    alpha = 0.4,
    inherit.aes = FALSE
  ) +
  # geom_path(aes(x = datetime, y = abs(bathy))) +
  geom_path() +
  scale_x_datetime(
    # date_breaks = "6 hours",
    # date_labels = "+ %Hh",
    breaks = seq.POSIXt(
      as.POSIXct("2000-01-01 06:00:00", tz = "UTC"),
      as.POSIXct("2000-01-02 06:00:00", tz = "UTC"),
      by = "6 hour"
    ),
    labels = c("-30h", "-24h", "-18h", "-12h", "-6h"),
    position = "top"
  ) +
  scale_y_reverse(limits = rev(range(data_plot$depth))) +
  labs(x = "Time", y = "Depth (m)") +
  # facet_nested(trip + id ~ phase,
  facet_nested(id ~ phase,
    scales = "free_x",
    # independent = "x"
  ) +
  theme(
    strip.placement = "outside",
    legend.position = "top",
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank(),
    panel.grid.major = element_line(color = "grey90"),
    axis.line.x = element_line(
      color = "black",
      arrow = arrow(length = unit(0.2, "lines"), type = "closed")
    ),
    axis.line.y = element_blank(),
    axis.ticks.y = element_blank(),
    axis.text.y = element_blank(),
    axis.title.y = element_blank(),
    axis.title.x = element_blank()
  )

fig_label_x <-
  ggplot(data.frame(l = "Time", x = 1, y = 1)) +
  geom_text(aes(x, y, label = l)) +
  theme_void() +
  coord_cartesian(clip = "off")

Code

guide_area() / fig_label_x / (departure_dp | arrival_dp) +
  plot_layout(
    guides = "collect",
    heights = c(0.5, 0.7, 6)
  ) &
  theme(
    legend.position = "top",
    plot.margin = unit(c(0, 0, 0, 0), "npc")
  )

Figure 1: Snapshot of departure and arrival dive profiles for six seals. Departure and arrival were constrained to 36 hours. Shading indicates night and coloration indicates dive type.

Notes

• TODO: Add three dots between Departure and Arrival

4.3 Figure 2

Let’s create a table specific for this figure that only contains for each individual:

• the first dive
• the last dive
• all benthic dives

Code

# let's add a column with the local_time
data_windrose <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # then by individual
  group_by(id) %>%
  # keep only the first, last date, but also all benthic dives
  filter(date == min(date) |
    date == max(date)) %>%
  # ungroup
  ungroup() %>%
  # change colnames to make getSunlightTimes works
  mutate(
    date = as.Date(date_tz),
    lat = Lat,
    lon = Long
  ) %>%
  # calculate hour of sunset and hour of sunrise
  mutate(getSunlightTimes(
    data = .,
    keep = c("sunrise", "sunset"),
    tz = tz_lookup_coords(
      lat = Lat,
      lon = Long,
      method = "accurate"
    )
  )) %>%
  # extract hours of sunset an sunrise
  mutate(
    h_sunrise = round(hour(sunrise) +
      minute(sunrise) / 60 +
      second(sunrise) / (60 * 60)),
    h_sunset = round(hour(sunset) +
      minute(sunset) / 60 +
      second(sunset) / (60 * 60))
  )

Code

# interesting circular time series
# https://campus.datacamp.com/courses/fraud-detection-in-r/introduction-motivation?ex=4

# for departure
windrose_departure_pb <- data_windrose %>%
  group_by(id) %>%
  filter(date_tz == min(date_tz) & trip == "Post Breeding") %>%
  mutate(time = as.numeric(format(date_tz, format = "%H"))) %>%
  group_by(time, trip) %>%
  summarise(nb_ind_hour = n(), .groups = "drop") %>%
  mutate(origin = "Departure") %>%
  windrose(.,
    grid = c(5, 10, 15),
    facet_wrap = "origin",
    hour_sunrise = data_windrose %>%
      group_by(id) %>%
      filter(date_tz == min(date_tz) &
        trip == "Post Breeding") %>%
      pull(h_sunrise) %>%
      median(.),
    hour_sunset = data_windrose %>%
      group_by(id) %>%
      filter(date_tz == min(date_tz) &
        trip == "Post Breeding") %>%
      pull(h_sunset) %>%
      median(.)
  ) +
  theme(
    axis.title = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  ) +
  # geom_vline(xintercept = 6,
  #            linetype = "dashed") +
  geom_text(
    data = data.frame(
      x = rep(8, 3),
      y = c(6, 11, 16),
      text = c(5, 10, 15)
    ),
    aes(x = x, y = y, label = text),
    size = 2
  )

windrose_departure_pm <- data_windrose %>%
  group_by(id) %>%
  filter(date_tz == min(date_tz) & trip == "Post Molting") %>%
  mutate(time = as.numeric(format(date_tz, format = "%H"))) %>%
  group_by(time, trip) %>%
  summarise(nb_ind_hour = n(), .groups = "drop") %>%
  mutate(origin = "Departure") %>%
  windrose(.,
    grid = c(5, 10, 15),
    hour_sunrise = data_windrose %>%
      group_by(id) %>%
      filter(date_tz == min(date_tz) &
        trip == "Post Molting") %>%
      pull(h_sunrise) %>%
      median(.),
    hour_sunset = data_windrose %>%
      group_by(id) %>%
      filter(date_tz == min(date_tz) &
        trip == "Post Molting") %>%
      pull(h_sunset) %>%
      median(.)
  ) +
  theme(
    axis.title = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  ) #+
# geom_vline(xintercept = 6,
#            linetype = "dashed") +
# geom_text(data = data.frame(
#   x = rep(8,3),
#   y = c(6, 11, 16),
#   text = c(5, 10, 15)
# ),
# aes(x = x, y = y, label = text),
# size = 2)

# for arrival
windrose_arrival_pb <- data_windrose %>%
  group_by(id) %>%
  filter(date_tz == max(date_tz) & trip == "Post Breeding") %>%
  mutate(time = as.numeric(format(date_tz, format = "%H"))) %>%
  group_by(time, trip) %>%
  summarise(nb_ind_hour = n(), .groups = "drop") %>%
  mutate(origin = "Arrival") %>%
  windrose(.,
    grid = c(5, 10, 15),
    legend_position = "top",
    facet_grid = T,
    hour_sunrise = data_windrose %>%
      group_by(id) %>%
      filter(date_tz == max(date_tz) &
        trip == "Post Breeding") %>%
      pull(h_sunrise) %>%
      median(.),
    hour_sunset = data_windrose %>%
      group_by(id) %>%
      filter(date_tz == max(date_tz) &
        trip == "Post Breeding") %>%
      pull(h_sunset) %>%
      median(.)
  ) +
  theme(
    axis.title = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  ) #+
#   geom_text(data = data.frame(
#   x = rep(8,3),
#   y = c(6, 11, 16),
#   text = c(5, 10, 15)
# ),
# aes(x = x, y = y, label = text),
# size = 2)

windrose_arrival_pm <- data_windrose %>%
  group_by(id) %>%
  filter(date_tz == max(date_tz) & trip == "Post Molting") %>%
  mutate(time = as.numeric(format(date_tz, format = "%H"))) %>%
  group_by(time, trip) %>%
  summarise(nb_ind_hour = n(), .groups = "drop") %>%
  mutate(origin = "Arrival") %>%
  windrose(.,
    grid = c(5, 10, 15),
    legend_position = "top",
    facet_wrap = "trip",
    hour_sunrise = data_windrose %>%
      group_by(id) %>%
      filter(date_tz == max(date_tz) &
        trip == "Post Molting") %>%
      pull(h_sunrise) %>%
      median(.),
    hour_sunset = data_windrose %>%
      group_by(id) %>%
      filter(date_tz == max(date_tz) &
        trip == "Post Molting") %>%
      pull(h_sunset) %>%
      median(.)
  ) +
  theme(
    axis.title = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  ) #+
#   geom_text(data = data.frame(
#   x = rep(8,3),
#   y = c(6, 11, 16),
#   text = c(5, 10, 15)
# ),
# aes(x = x, y = y, label = text),
# size = 2)

label_x <- ggplot(data.frame(l = "Time (hours)", x = 1, y = 1)) +
  geom_text(aes(x, y, label = l)) +
  theme_void() +
  coord_cartesian(clip = "off")

# combine all
(windrose_departure_pb + windrose_arrival_pb) /
  (windrose_departure_pm + windrose_arrival_pm) /
  label_x +
  theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
  plot_layout(
    heights = c(2, 2, 0),
    guides = "collect"
  ) &
  # plot_annotation(tag_levels = c("A")) &
  # theme(plot.tag.position = c(-0.03, 0.8),
  #       plot.tag = element_text(size = 8, hjust = 0, vjust = 0)) +
  theme(legend.position = "top")

Figure 2: Circular histogram plots display the local timing (in hours) for the first and last recorded dives (respectively departure and arrival) of adult female northern elephant seals (n = 403) during their Post Breeding and Post Molting foraging trips.

Note

• n = 401, because we removed individuals without location data (filter(!is.na(Lat)))

Code

# departure PM
pm_dep_mean <- psych::circadian.mean(
  data_windrose %>%
    group_by(id) %>%
    filter(date_tz == min(date_tz) & trip == "Post Molting") %>%
    mutate(time = hour(date_tz) +
      minute(date_tz) / 60 +
      second(date_tz) / (60 * 60)) %>% pull(time)
)

# departure PB
pb_dep_mean <- psych::circadian.mean(
  data_windrose %>%
    group_by(id) %>%
    filter(date_tz == min(date_tz) & trip == "Post Breeding") %>%
    mutate(time = hour(date_tz) +
      minute(date_tz) / 60 +
      second(date_tz) / (60 * 60)) %>% pull(time)
)

# arrival PM
pm_arr_mean <- psych::circadian.mean(
  data_windrose %>%
    group_by(id) %>%
    filter(date_tz == max(date_tz) & trip == "Post Molting") %>%
    mutate(time = hour(date_tz) +
      minute(date_tz) / 60 +
      second(date_tz) / (60 * 60)) %>% pull(time)
)

# arrival PB
pb_arr_mean <- psych::circadian.mean(
  data_windrose %>%
    group_by(id) %>%
    filter(date_tz == max(date_tz) & trip == "Post Breeding") %>%
    mutate(time = hour(date_tz) +
      minute(date_tz) / 60 +
      second(date_tz) / (60 * 60)) %>% pull(time)
)

Code

# functio to convert time in radian (probably already exist in some package...)
hour_time_2_radian <- function(x, high = 24, low = 0) {
  # https://pingouin-stats.org/build/html/generated/pingouin.convert_angles.html#pingouin.convert_angles
  ptp <- high - low
  rad <- x * (2 * pi) / ptp
  rad <- (rad + pi) %% (2 * pi) - pi
  return(rad)
}

Code

# departure
test_departure_pm <- CircStats::r.test(
  data_windrose %>%
    group_by(id) %>%
    filter(date_tz == min(date_tz) & trip == "Post Molting") %>%
    mutate(time = hour(date_tz) +
      minute(date_tz) / 60 +
      second(date_tz) / (60 * 60)) %>%
    pull(time) %>%
    hour_time_2_radian(.)
)
test_departure_pb <- CircStats::r.test(
  data_windrose %>%
    group_by(id) %>%
    filter(date_tz == min(date_tz) & trip == "Post Breeding") %>%
    mutate(time = hour(date_tz) +
      minute(date_tz) / 60 +
      second(date_tz) / (60 * 60)) %>%
    pull(time) %>%
    hour_time_2_radian(.)
)

# arrival
test_arrival_pm <- CircStats::r.test(
  data_windrose %>%
    group_by(id) %>%
    filter(date_tz == max(date_tz) & trip == "Post Molting") %>%
    mutate(time = hour(date_tz) +
      minute(date_tz) / 60 +
      second(date_tz) / (60 * 60)) %>%
    pull(time) %>%
    hour_time_2_radian(.)
)
test_arrival_pb <- CircStats::r.test(
  data_windrose %>%
    group_by(id) %>%
    filter(date_tz == max(date_tz) & trip == "Post Breeding") %>%
    mutate(time = hour(date_tz) +
      minute(date_tz) / 60 +
      second(date_tz) / (60 * 60)) %>%
    pull(time) %>%
    hour_time_2_radian(.)
)

Code

data.frame(
  Departure = c(
    paste0(
      round(pb_dep_mean, 1),
      " (R = ",
      round(test_departure_pb$r.bar, 2),
      ", p-value ",
      ifelse(
        round(test_departure_pb$p.value, 2) == 0,
        "< 0.001",
        paste("= ", round(test_departure_pb$p.value, 2))
      ),
      ")"
    ),
    paste0(
      round(pm_dep_mean, 1),
      " (R = ",
      round(test_departure_pm$r.bar, 2),
      ", p-value ",
      ifelse(
        round(test_departure_pm$p.value, 2) == 0,
        "< 0.001",
        paste("= ", round(test_departure_pm$p.value, 2))
      ),
      ")"
    )
  ),
  Arrival = c(
    paste0(
      round(pb_arr_mean, 1),
      " (R = ",
      round(test_arrival_pb$r.bar, 2),
      ", p-value ",
      ifelse(
        round(test_arrival_pb$p.value, 2) == 0,
        "< 0.001",
        paste("= ", round(test_arrival_pb$p.value, 2))
      ),
      ")"
    ),
    paste0(
      round(pm_arr_mean, 1),
      " (R = ",
      round(test_arrival_pm$r.bar, 2),
      ", p-value ",
      ifelse(
        round(test_arrival_pm$p.value, 2) == 0,
        "< 0.001",
        paste("= ", round(test_arrival_pm$p.value, 2))
      ),
      ")"
    )
  ),
  rownames = c("Post Breeding", "Post Molting")
) %>%
  gt(rowname_col = "rownames") %>%
  cols_label(
    Arrival = "Hour of Arrival",
    Departure = "Hour of Departure"
  )

Table 3:

(Not meant to be in the paper) Average hour of departure and arrival for post breeding and post molting trip with the result of the associated Rayleigh test.

	Hour of Departure	Hour of Arrival
Post Breeding	19 (R = 0.2, p-value < 0.001)	15.8 (R = 0.04, p-value = 0.65)
Post Molting	19.2 (R = 0.22, p-value < 0.001)	13.8 (R = 0.13, p-value = 0.06)

Good ressource to interpret Rayleigh test: https://stats.stackexchange.com/questions/433396/interpretation-of-rayleigh-test-results

4.4 Figure 3

Code

# with or without bathymetry
with_bathy <- T

# if bathy
if (with_bathy) {
  # check if background_ggoceanmaps exist
  if (file.exists("../export/background_ggoceanmap_zoom_in.rds")) {
    trip_zoom_in <- readRDS("../export/background_ggoceanmap_zoom_in.rds")
  } else {
    # using ggOceanMaps
    trip_zoom_in <- basemap(
      # limits = c(170, -110, 30, 59),
      limits = c(-155, -110, 30, 60),
      bathymetry = TRUE,
      shapefiles = "Arctic",
      rotate = TRUE,
      grid.col = NA
    )

    # Make the graticules:
    lims <- attributes(trip_zoom_in)$limits
    graticule <- sf::st_graticule(
      c(lims[1], lims[3], lims[2], lims[4]),
      crs = attributes(trip_zoom_in)$proj,
      lon = seq(-180, 180, 45),
      lat = seq(-90, 90, 10)
    )

    # Plot
    trip_zoom_in <- trip_zoom_in +
      geom_sf(data = graticule, color = "grey50")

    # reorder layers
    trip_zoom_in$layers <- trip_zoom_in$layers[c(1, 3, 2)]

    # save result
    saveRDS(trip_zoom_in, "../export/background_ggoceanmap_zoom_in.rds")
  }
} else {
  # using ggOceanMaps
  trip_zoom_in <- basemap(
    # limits = c(170, -110, 20, 59),
    limits = c(-155, -110, 30, 60),
    land.col = "black",
    # land.col = '#cdaa6d',
    bathymetry = FALSE,
    shapefiles = "Arctic",
    rotate = TRUE,
    grid.col = NA
  )

  # Make the graticules:
  lims <- attributes(trip_zoom_in)$limits
  graticule <- sf::st_graticule(
    c(lims[1], lims[3], lims[2], lims[4]),
    crs = attributes(trip_zoom_in)$proj,
    lon = seq(-180, 180, 45),
    lat = seq(-90, 90, 10)
  )

  # Plot
  trip_zoom_in <- trip_zoom_in +
    geom_sf(data = graticule, color = "grey70", linewidth = 0.3) +
    theme(
      panel.background = element_rect(fill = "#4292c6"),
      panel.ontop = FALSE
    )

  # reorder layers
  trip_zoom_in$layers <- trip_zoom_in$layers[c(2, 1)]
}

Code

# get predator distribution *.kmz (extraction using google earth)
# for calculation
Shark_area_coord <- read_sf("../export/shark.kmz") %>%
  smooth(., method = "ksmooth", smoothness = 0.5) %>%
  # st_crop(., xmin = 170, xmax = -110, ymin = 20, ymax = 59) %>%
  st_coordinates() %>%
  as_tibble() %>%
  select(lon = X, lat = Y)
Orca_area_coord <- read_sf("../export/orca.kmz") %>%
  smooth(., method = "ksmooth", smoothness = 0.5) %>%
  # st_crop(., xmin = 170, xmax = -110, ymin = 20, ymax = 59) %>%
  st_coordinates() %>%
  as_tibble() %>%
  select(lon = X, lat = Y)

# for plots
Overlap_area_coord_to_display <- st_intersection(
  read_sf("../export/shark.kmz") %>%
    smooth(., method = "ksmooth", smoothness = 0.5),
  read_sf("../export/orca.kmz") %>%
    smooth(., method = "ksmooth", smoothness = 0.5)
) %>%
  st_coordinates() %>%
  as_tibble() %>%
  select(lon = X, lat = Y)
Shark_area_coord_to_display <- st_difference(
  read_sf("../export/shark.kmz") %>%
    smooth(., method = "ksmooth", smoothness = 0.5),
  read_sf("../export/orca.kmz") %>%
    smooth(., method = "ksmooth", smoothness = 0.5)
) %>%
  # st_crop(., xmin = 170, xmax = -110, ymin = 20, ymax = 59) %>%
  st_coordinates() %>%
  as_tibble() %>%
  select(lon = X, lat = Y)
Orca_area_coord_to_display <- st_difference(
  read_sf("../export/orca.kmz") %>%
    smooth(., method = "ksmooth", smoothness = 0.5),
  read_sf("../export/shark.kmz") %>%
    smooth(., method = "ksmooth", smoothness = 0.5)
) %>%
  # st_crop(., xmin = 170, xmax = -110, ymin = 20, ymax = 59) %>%
  st_coordinates() %>%
  as_tibble() %>%
  select(lon = X, lat = Y)

Code

# add inside/outside each area
data_dive <- data_dive %>%
  mutate(
    shark_area = point.in.polygon(
      point.x = Long,
      point.y = Lat,
      pol.x = Shark_area_coord$lon,
      pol.y = Shark_area_coord$lat
    ),
    orca_area = point.in.polygon(
      point.x = Long,
      point.y = Lat,
      pol.x = Orca_area_coord$lon,
      pol.y = Orca_area_coord$lat
    ),
    overlap_area = as.numeric(if_else(shark_area == 1 &
      orca_area == 1, T, F)),
    total_area = as.numeric(if_else(shark_area == 0 &
      orca_area == 0, F, T))
  )

Code

# data_dive_longer <- pivot_longer(
#   data_dive %>%
#     # keep data with locations
#     filter(!is.na(Lat)),
#   c(
#     "shark_area",
#     "orca_area",
#     "overlap_area",
#     "total_area"
#   ),
#   names_to = "areas",
#   values_to = "in_out"
# )
#
# data_descriptive_ind <- data_dive_longer %>%
#   filter(DiveTypeName == "Benthic") %>%
#   group_by(areas, id) %>%
#   summarise(
#     n_tot = n(),
#     n_inside = length(id[in_out == 1]),
#     n_outside = length(id[in_out == 0]),
#     .groups = "drop"
#   ) %>%
#   mutate(
#     perc_inside = n_inside / n_tot,
#     perc_outside = n_outside / n_tot
#   )
#
# data_descriptive_area <- data_descriptive_ind %>%
#   group_by(areas) %>%
#   summarise(
#     mean_perc_inside = wtd.mean(perc_inside, n_tot),
#     mean_perc_outside = wtd.mean(perc_outside, n_tot),
#     sd_perc_inside = sqrt(wtd.var(perc_inside, n_tot)),
#     sd_perc_outside = sqrt(wtd.var(perc_outside, n_tot)),
#     # based from here: https://stats.stackexchange.com/questions/25895/computing-standard-error-in-weighted-mean-estimation
#     se_perc_inside = sqrt(wtd.var(perc_inside, n_tot)*sum((n_tot/sum(n_tot))^2)),
#     se_perc_outside = sqrt(wtd.var(perc_outside, n_tot)*sum((n_tot/sum(n_tot))^2)),
#   ) %>%
#   pivot_longer(
#     cols = ends_with("inside") |
#       ends_with("outside"),
#     names_to = c("stat", "in_out"),
#     names_pattern = "(.*)_perc_(.*)"
#   ) %>%
#   pivot_wider(
#     id_cols = c("areas", "in_out"),
#     names_from = "stat",
#     values_from = "value"
#   ) %>%
#   mutate(areas = factor(areas,
#                         level = c(
#                           "overlap_area",
#                           "orca_area",
#                           "shark_area",
#                           "total_area"
#                         )
#   ))
#
# data_stat_area <- data_dive_longer %>%
#   filter(DiveTypeName == "Benthic") %>%
#   group_by(areas, in_out) %>%
#   summarise(nb_area = n(), .groups = "drop_last") %>%
#   mutate(nb_total = sum(nb_area)) %>%
#   summarise(rstatix::prop_test(x = nb_area, n = nb_total))

Code

# inside_bar_plot = data_descriptive_area %>%
#   mutate(
#     mean = if_else(in_out == "outside", -mean, mean),
#     sd = if_else(in_out == "outside", -sd, sd),
#     se = if_else(in_out == "outside", -se, se)
#   ) %>%
#   filter(in_out == "inside") %>%
#   ggplot(aes(x = mean, y = areas, fill = areas)) +
#   facet_grid(in_out ~ ., scales = "free_y") +
#   geom_col(show.legend = F) +
#   scale_fill_manual(values = c("#9A6BEC", "#E66100", "#E1E826", "#004D40")) +
#   new_scale_fill() +
#   geom_col(
#     show.legend = F,
#     mapping = aes(
#       x = mean,
#       y = areas,
#       fill = in_out),
#     alpha = 0.3
#   ) +
#   scale_fill_manual(values = c("grey70")) +
#   geom_errorbar(mapping = aes(xmin = -se,
#                               xmax = se),
#                 # reduce size of horizontal bar of the errorbar
#                 width = 0.05) +
#   coord_flip(xlim = c(0, max(data_descriptive_area$mean))) +
#   scale_x_continuous(
#     # breaks = c(-0.6,-0.4,-0.2, 0.2, 0.4, 0.6),
#     expand = expansion(mult = c(0,0.1), add = 0),
#     labels = function(x) {
#       scales::percent(abs(x), 1)
#     }
#   ) +
#   scale_y_discrete(
#     labels = function(x) {
#       lapply(x, function(x)
#         (str_split(str_to_title(x), "_"))[[1]][1])
#     }
#   ) +
#   theme(
#     panel.spacing.y = unit(0, "mm"),
#     panel.grid.major = element_blank(),
#     panel.grid.minor = element_blank(),
#     panel.background = element_blank(),
#     axis.line.y = element_line(arrow = arrow(
#       angle = 20,
#       length = unit(.10, "inches"),
#       type = "closed"
#     )),
#     axis.title = element_blank(),
#     axis.ticks.length.x = unit(0, units = "mm"),
#     axis.text.x = element_blank()
#   )
#
# outside_bar_plot = data_descriptive_area %>%
#   mutate(
#     mean = if_else(in_out == "outside",-mean, mean),
#     sd = if_else(in_out == "outside",-sd, sd),
#     se = if_else(in_out == "outside",-se, se)
#   ) %>%
#   filter(in_out == "outside") %>%
#   ggplot(aes(x = mean, y = areas, fill = areas)) +
#   facet_grid(in_out ~ ., scales = "free_y") +
#   geom_col(show.legend = F) +
#   scale_fill_manual(values = c("#9A6BEC", "#E66100", "#E1E826", "#004D40")) +
#   new_scale_fill() +
#   geom_col(show.legend = F, mapping = aes(x = mean,
#                                           y = areas,
#                                           fill = in_out),
#                                           alpha = 0.3) +
#   scale_fill_manual(values = c("grey30")) +
#   geom_errorbar(
#     mapping = aes(
#       xmin = -se,
#       xmax = se
#     ),
#     # reduce size of horizontal bar of the errorbar
#     width = 0.05
#   ) +
#   coord_flip(xlim = c(-max(data_descriptive_area$mean), 0)) +
#   scale_x_continuous(
#     # breaks = c(-0.6, -0.4, -0.2, 0.2, 0.4, 0.6),
#     expand = expansion(mult = c(0.1,0), add = 0),
#     labels = function(x) {
#       scales::percent(abs(x), 1)
#     }
#   ) +
#   scale_y_discrete(
#     labels = function(x) {
#       lapply(x, function(x) (str_split(str_to_title(x), "_"))[[1]][1])
#     }
#   ) +
#   labs(y = "Areas") +
#   theme(panel.spacing.y = unit(0, "mm"),
#         panel.grid.major = element_blank(),
#         panel.grid.minor = element_blank(),
#         panel.background = element_blank(),
#         axis.line.y = element_line(arrow = arrow(
#           angle = 20,
#           length = unit(.10, "inches"),
#           type = "closed",
#           ends = "first"
#         )),
#         axis.ticks.length.x = unit(0,units = "mm"),
#         axis.title.y = element_blank())
#
# label_x_bar_plot = ggplot(data.frame(l = "Proportion of benthic dives",
#                                      x = 1,
#                                      y = 1)) +
#   geom_text(aes(x, y, label = l, angle = 90)) +
#   theme(aspect.ratio = 10,
#         axis.title = element_blank(),
#         axis.text = element_blank(),
#         axis.ticks = element_blank()) +
#   theme_void()

Code

# color palette
# https://davidmathlogic.com/colorblind/#%23D81B60-%231E88E5-%23FFC107-%23004D40
# combine areas
areas_coord <- rbind(
  Orca_area_coord_to_display %>% mutate(Predator = "orca"),
  Shark_area_coord_to_display %>% mutate(Predator = "shark"),
  Overlap_area_coord_to_display %>% mutate(Predator = "overlap")
) %>%
  mutate(Predator = factor(Predator, levels = c("orca", "shark", "overlap")))

# data use to compute kernel density estimation
df_kernel_dens_benthic <- data_dive %>%
  # only with location data
  filter(!is.na(Lat) & DiveTypeName == "Benthic") %>%
  # select only the required columns
  select(lat = Lat, long = Long, id, DiveTypeName) %>%
  # create col to nicely display dive type
  mutate(origin = paste(DiveTypeName, "dives"))

# transform data into sf object
df_kernel_dens_sf_benthic <- st_as_sf(
  df_kernel_dens_benthic,
  coords = c("long", "lat"),
  crs = st_crs(4326)
)

# make it's group_by origin
df_kernel_dens_sf_benthic <- group_by(df_kernel_dens_sf_benthic, origin)

# kernel density estimation
df_kernel_dens_sf_kde_benthic <- st_kde(df_kernel_dens_sf_benthic,
  H = diag(c(
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf_benthic)[, 1]
    ),
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf_benthic)[, 2]
    )
  ) / 4)^2
)

# https://loading.io/color/feature/BuGn-9/
cols <- c(
  "99%" = "#ccece6",
  "95%" = "#99d8c9",
  "80%" = "#41ae76",
  "50%" = "#00441b"
)

p1 <- trip_zoom_in +
  # remove legend
  guides(fill = "none") +
  # new scale
  new_scale_fill() +
  # kernel
  geom_sf(
    data = st_get_contour( # geospatial kernel
      df_kernel_dens_sf_kde_benthic,
      # probabilities
      cont = c(50, 80, 95, 99)
    ),
    # display
    aes(fill = label_percent(contlabel)),
    alpha = 0.7
  ) +
  # same colour bar
  # scale_fill_viridis_d(option = "magma", direction = -1) +
  # scale_fill_brewer(palette = "heat")
  # scale_fill_brewer(palette = "BuGn") +
  scale_fill_manual(values = cols) +
  # scale_fill_grey() +
  # legend
  labs(fill = "Probs") +
  # legend at the top
  theme(legend.position = "top") +
  # set limit from predator_map$layers[[1]]$data
  coord_sf(
    xlim = c(
      st_bbox(trip_zoom_in$layers[[1]]$data)$xmin,
      st_bbox(trip_zoom_in$layers[[1]]$data)$xmax
    ),
    ylim = c(
      st_bbox(trip_zoom_in$layers[[1]]$data)$ymin,
      st_bbox(trip_zoom_in$layers[[1]]$data)$ymax
    )
  )


# Extract the size legend - leg1
leg1 <- gtable_filter(ggplot_gtable(ggplot_build(p1)), "guide-box")

p2 <- trip_zoom_in +
  # tweak legend
  guides(fill = guide_legend(
    override.aes = list(
      colour = "grey60"
    )
  )) +
  theme(
    legend.key = element_rect(color = "grey60"),
    legend.text = element_text(colour = "grey20"),
    legend.title = element_text(colour = "grey20"),
    legend.justification = c("right", "top"),
    legend.position = c(0.99, 0.99),
    legend.background = element_rect(fill = alpha("white", 0))
  ) +
  labs(x = "Longitude", y = "Latitude") +
  # new scale
  new_scale_fill() +
  # kernel
  geom_sf(
    data = st_get_contour( # geospatial kernel
      df_kernel_dens_sf_kde_benthic,
      # probabilities
      cont = c(50, 80, 95, 99)
    ),
    # display
    aes(fill = label_percent(contlabel)),
    alpha = 0.7
  ) +
  # same colour bar
  # scale_fill_viridis_d(option = "magma", direction = -1) +
  # scale_fill_brewer(palette = "heat")
  # scale_fill_brewer(palette = "BuGn") +
  scale_fill_manual(values = cols) +
  # scale_fill_grey() +
  # remove legend
  guides(fill = "none") +
  # # no display alpha
  # guides(alpha = "none") +
  # new scale
  new_scale_fill() +
  # add spatial polygon
  geom_spatial_polygon(
    data = areas_coord,
    aes(
      x = lon,
      y = lat,
      fill = Predator
    ),
    alpha = 0.35,
    crs = 4326
  ) +
  # color palette
  scale_fill_manual(values = c("darkred", "#E1E826", "#E66100")) +
  scale_color_manual(values = c("darkred", "#E1E826", "#E66100")) +
  # remove legend for areas
  guides(fill = "none", col = "none") +
  # set limit from predator_map$layers[[1]]$data
  coord_sf(
    xlim = c(
      st_bbox(trip_zoom_in$layers[[1]]$data)$xmin,
      st_bbox(trip_zoom_in$layers[[1]]$data)$xmax
    ),
    ylim = c(
      st_bbox(trip_zoom_in$layers[[1]]$data)$ymin,
      st_bbox(trip_zoom_in$layers[[1]]$data)$ymax
    )
  )


# reordering layers (continent on top)
p2$layers <- if (with_bathy) {
  p2$layers[c(1, 2, 5, 4, 3)]
} else {
  p2$layers[c(1, 4, 3, 2)]
}

# convert into ggplot
predator_map <- as_ggplot(leg1) / p2 +
  plot_layout(
    ncol = 1,
    heights = c(1, 12)
  )

Code

# load image icon
shark_img <- readPNG("../img/shark_custom.png", native = T)
orca_img <- readPNG("../img/orca_custom.png", native = T)

# draw the map
predator_map +
  inset_element(
    p = orca_img,
    left = 0.15,
    bottom = 0.82,
    right = 0.25,
    top = 0.92
  ) +
  theme(
    rect = element_rect(fill = "transparent"),
    plot.background = element_rect(color = NA)
  ) +
  inset_element(
    p = shark_img,
    left = 0.15,
    bottom = 0.3,
    right = 0.25,
    top = 0.4
  ) +
  theme(
    rect = element_rect(fill = "transparent"),
    plot.background = element_rect(color = NA)
  )

Figure 3: Spatial overlap of benthic dives performed by northern elephant seals (n = 401), and the shark (yellow polygon) and orca (red polygon) distribution area along the west coast (adapted from Jorgenson, et al. 2019) with a kernel density estimation representing the occurrence of all the benthic dives.

Code

# # layout
# layout <- "
# #AAAAAAAAAAAAAAAAAAAAAAA
# #AAAAAAAAAAAAAAAAAAAAAAA
# #AAAAAAAAAAAAAAAAAAAAAAA
# #AAAAAAAAAAAAAAAAAAAAAAA
# #AAAAAAAAAAAAAAAAAAAAAAA
# BCCCCCCCCCCCCCCCCCCCCCCC
# BDDDDDDDDDDDDDDDDDDDDDDD
# "
#
# # predator_map$layers[[length(predator_map$layers)]]$aes_params$fill_new_new = "white" #cdaa6d
# # predator_map$layers[[length(predator_map$layers)]]$aes_params$colour = "black"
#
# # display
# (predator_map + theme(plot.margin = unit(c(0, 0, 0, 0), "pt"))) +
#   (label_x_bar_plot + theme(plot.margin = unit(c(0, 0, 0, 0), "pt"))) +
#   (inside_bar_plot + theme(plot.margin = unit(c(0, 0, 0, 0), "pt"))) +
#   (outside_bar_plot + theme(plot.margin = unit(c(0, 0, 0, 0), "pt"))) +
#   plot_layout(design = layout) #+
#   #plot_annotation(tag_levels = list(c("A","","B","")))

Notes

• ?
• n = 401, because we removed individuals without location data and those that do not have any benthic dives (filter(!is.na(Lat) & DiveTypeName == "Benthic"))

4.5 Figure 4

Code

# https://mycolor.space/?hex=%23008C49&sub=1
# build palette color
cols_ind_1 <- c("#008c49", "#00c9d3", "#bffcfd", "#458081")
# test for color blind
# => https://davidmathlogic.com/colorblind/#%23008C49-%2300C9D3-%23BFFCFD-%23458081
cols_ind_2 <- c("#008c49", "#002e00", "#00c9d3", "#deb887")
# test for color blind
# => https://davidmathlogic.com/colorblind/#%23008C49-%23002E00-%2300C9D3-%23DEB887

Code

# main plot
bathy_prop <- data_dive %>%
  # keep only rows with location data
  filter(!is.na(Lat)) %>%
  # keep only negative bathy
  filter(bathy < 0) %>%
  # and remove outliers
  filter(bathy > -6000) %>%
  # positive bathy
  mutate(bathy = abs(bathy)) %>%
  # create class of bathymetry
  mutate(bathy_class = cut(
    bathy,
    seq(0, 6000, 400),
    ordered_result = T,
    dig.lab = 4
  )) %>%
  # calculate by bath_class and animal
  group_by(bathy_class, DiveTypeName) %>%
  # the number of dives
  summarise(N = n(), .groups = "drop") %>%
  # and divide by the total number of dives per bathy_class
  group_by(bathy_class) %>%
  # to get the Proportion of different dive types per bathy_class
  mutate(Proportion = N / sum(N)) %>%
  # ungroup => not required but that let the dataset clean
  ungroup() %>%
  # then plot
  ggplot(aes(x = bathy_class, y = Proportion)) +
  # the area
  geom_area(aes(fill = DiveTypeName, group = DiveTypeName)) +
  # orientation of x labels
  guides(x = guide_axis(angle = 45)) +
  # format y axis
  scale_y_continuous(
    labels = function(x) {
      paste0(x * 100, "%")
    }
  ) +
  labs(
    x = "Bathymetry  (m)",
    y = "Dive type proportion (%)"
  ) +
  # scale_fill_viridis(
  #   "Dive Type",
  #   option = "plasma",
  #   discrete = T,
  #   direction = -1
  # ) +
  scale_fill_manual("Dive Type",
    values = cols_ind_1
  ) +
  # position legend at the top
  theme(
    legend.position = "top",
    axis.title.x = element_text(margin = margin(t = 15)),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.line = element_line(arrow = arrow(
      angle = 20,
      length = unit(.10, "inches"),
      type = "closed"
    )),
    plot.tag.position = c(-0.03, 1)
  )

# # second plot
# bathy_hist <- data_dive %>%
#   # keep only negative bathy
#   filter(bathy < 0) %>%
#   # and remove outliers
#   filter(bathy > -6000) %>%
#   # create class of bathymetry
#   mutate(bathy_class = fct_rev(cut(
#     bathy,
#     seq(-6000, 0, 400),
#     ordered_result = T,
#     dig.lab = 4
#   ))) %>%
#   # calculate by bath_class and animal
#   group_by(bathy_class, DiveTypeName) %>%
#   # the number of dives
#   summarise(N = n(), .groups = "drop") %>%
#   # ggplot
#   ggplot(aes(
#     x = bathy_class,
#     y = N,
#     fill = DiveTypeName,
#     width = 0.5
#   )) +
#   geom_bar(
#     stat = "identity",
#     position = "dodge"
#   ) +
#   scale_fill_viridis(
#     option = "plasma",
#     discrete = T,
#     direction = -1
#   ) +
#   theme_void() +
#   theme(legend.position = "top")
#
# # the actual plot
# bath_plot <- bathy_hist / bathy_prop + plot_layout(heights = c(1, 10))

Code

# main plot
dist_coast_prop <- data_dive %>%
  # filter out animal without dist_coast
  filter(dist_coast > 0) %>%
  # remove outliers
  filter(dist_coast * 1000 * 111 <= 1900) %>%
  # create class of bathymetry
  mutate(dist_class = cut(
    # convert decimal degree/1000
    dist_coast * 1000 * 111,
    seq(0, 1900, 100),
    ordered_result = T,
    dig.lab = 4
  )) %>%
  # calculate by bath_class and animal
  group_by(dist_class, DiveTypeName) %>%
  # the number of dives
  summarise(N = n(), .groups = "drop") %>%
  # and divide by the total number of dives per dist_class
  group_by(dist_class) %>%
  # to get the Proportion of different dive types per dist_class
  mutate(Proportion = N / sum(N)) %>%
  # ungroup => not required but that let the dataset clean
  ungroup() %>%
  # then plot
  ggplot(aes(x = dist_class, y = Proportion)) +
  # the area
  geom_area(aes(fill = DiveTypeName, group = DiveTypeName)) +
  # orientation of x labels
  guides(x = guide_axis(angle = 45)) +
  # format y axis
  scale_y_continuous(
    labels = function(x) {
      paste0(x * 100, "%")
    }
  ) +
  labs(
    x = "Distance from the coast (km)",
    y = "Dive type proportion"
  ) +
  # scale_fill_viridis(
  #   "Dive Type",
  #   option = "plasma",
  #   discrete = T,
  #   direction = -1
  # ) +
  scale_fill_manual("Dive Type",
    values = cols_ind_1
  ) +
  # position legend at the top
  theme(
    legend.position = "none",
    axis.title.x = element_text(margin = margin(t = 15)),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.line = element_line(arrow = arrow(
      angle = 20,
      length = unit(.10, "inches"),
      type = "closed"
    )),
    plot.tag.position = c(-0.03, 1)
  )

# # the second plot
# dist_coast_hist <- data_dive %>%
#   # filter out animal without dist_coast
#   filter(dist_coast > 0) %>%
#   # remove outliers
#   filter(dist_coast * 1000 * 111 <= 1900) %>%
#   # create class of bathymetry
#   mutate(dist_class = cut(
#     # convert decimal degree/1000
#     dist_coast * 1000 * 111,
#     seq(0, 2000, 100),
#     ordered_result = T,
#     dig.lab = 4
#   )) %>%
#   # calculate by bath_class and animal
#   group_by(dist_class, DiveTypeName) %>%
#   # the number of dives
#   summarise(N = n(), .groups = "drop") %>%
#   # ggplot
#   ggplot(aes(
#     x = dist_class,
#     y = N,
#     fill = DiveTypeName,
#     width = 0.5
#   )) +
#   geom_bar(
#     stat = "identity",
#     position = "dodge"
#   ) +
#   scale_fill_viridis(
#     option = "plasma",
#     discrete = T,
#     direction = -1
#   ) +
#   theme_void() +
#   theme(legend.position = "top")
#
# the actual plot
# dist_coast_plot <- dist_coast_hist / dist_coast_prop + plot_layout(heights = c(1, 10))

Code

# the main plot
dist_dep_prop <- data_dive %>%
  # filter out animal without dist_coast
  filter(dist_dep > 0) %>%
  # remove outliers
  filter(dist_dep / 1000 < 4200) %>%
  # create class of bathymetry
  mutate(dist_class = cut(
    dist_dep / 1000,
    seq(0, 4200, 300),
    ordered_result = T,
    dig.lab = 4
  )) %>%
  # calculate by bath_class and animal
  group_by(dist_class, DiveTypeName) %>%
  # the number of dives
  summarise(N = n(), .groups = "drop") %>%
  # and divide by the total number of dives per dist_class
  group_by(dist_class) %>%
  # to get the Proportion of different dive types per dist_class
  mutate(Proportion = N / sum(N)) %>%
  # ungroup => not required but that let the dataset clean
  ungroup() %>%
  # then plot
  ggplot(aes(x = dist_class, y = Proportion)) +
  # the area
  geom_area(aes(fill = DiveTypeName, group = DiveTypeName)) +
  # orientation of x labels
  guides(x = guide_axis(angle = 45)) +
  # format y axis
  scale_y_continuous(
    labels = function(x) {
      paste0(x * 100, "%")
    }
  ) +
  labs(
    x = "Distance from departure (km)",
    y = "Dive type proportion"
  ) +
  # scale_fill_viridis(
  #   "Dive Type",
  #   option = "plasma",
  #   discrete = T,
  #   direction = -1
  # ) +
  scale_fill_manual("Dive Type",
    values = cols_ind_1
  ) +
  # position legend at the top
  theme(
    legend.position = "none",
    axis.title.x = element_text(margin = margin(t = 15)),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.line = element_line(arrow = arrow(
      angle = 20,
      length = unit(.10, "inches"),
      type = "closed"
    )),
    plot.tag.position = c(-0.03, 1)
  )

# # the second plot
# dist_dep_hist <- data_dive %>%
#   # filter out animal without dist_coast
#   filter(dist_coast > 0) %>%
#   # remove outliers
#   filter(dist_dep / 1000 < 4200) %>%
#   # create class of bathymetry
#   mutate(dist_class = cut(
#     dist_dep / 1000,
#     seq(0, 4200, 300),
#     ordered_result = T,
#     dig.lab = 4
#   )) %>%
#   # calculate by bath_class and animal
#   group_by(dist_class, DiveTypeName) %>%
#   # the number of dives
#   summarise(N = n(), .groups = "drop") %>%
#   # ggplot
#   ggplot(aes(
#     x = dist_class,
#     y = N,
#     fill = DiveTypeName,
#     width = 0.5
#   )) +
#   geom_bar(
#     stat = "identity",
#     position = "dodge"
#   ) +
#   scale_fill_viridis(
#     option = "plasma",
#     discrete = T,
#     direction = -1
#   ) +
#   theme_void() +
#   theme(legend.position = "top")
#
# # the actual plot
# dist_dep_plot <- dist_dep_hist / dist_dep_prop + plot_layout(heights = c(1, 10))

Code

# y axis label position
ylab <- bathy_prop$labels$y
bathy_prop$labels$y <- dist_coast_prop$labels$y <- dist_dep_prop$labels$y <- " "

Code

# display
(guide_area() + bathy_prop + dist_coast_prop + dist_dep_prop) +
  theme(plot.margin = unit(c(0, .02, -1, .05), units = "npc")) +
  plot_annotation(tag_levels = "A") +
  plot_layout(nrow = 5, guides = "collect")
grid::grid.draw(grid::textGrob(ylab, x = 0.07, y = 0.55, rot = 90))

Figure 4: Proportion of each dive type performed by northern elephant seals in relation to (A) bathymetry, (B) distance from the coast, and (C) distance from departure (defined as the first recorded location).

Notes

• TODO: change the color code to match the color we attributed to benthic dives (greens?)
• n = 401, because we removed individuals without location data (filter(!is.na(Lat)))

4.6 Figure 5

Code

# the last plot
prop_at_sea_dt <- data_dive %>%
  group_by(id) %>%
  # time difference between the first date and the others (telling R to take diff b/w first date and all other dates)
  mutate(nb_days_departure = trunc(as.numeric(difftime(
    date, first(date),
    units = "days"
  )))) %>%
  # group by day/dive_type/seal
  group_by(nb_days_departure, DiveTypeName, id) %>%
  # count number of dives
  summarise(nb_daily_dives_type = n(), .groups = "drop") %>%
  # group by day/seal
  group_by(nb_days_departure, id) %>%
  # count the total number of dives per day/seal
  mutate(nb_daily_dives = sum(nb_daily_dives_type)) %>%
  # order by seals/date
  arrange(id, nb_days_departure) %>%
  # perform our calculation per seal
  group_by(id) %>%
  # calculate of the Proportion of time since departure
  mutate(perc_time_at_sea = round(
    nb_days_departure / max(nb_days_departure) * 100, 1
  )) %>%
  # filter on benthic dives
  filter(DiveTypeName == "Benthic") %>%
  # add class of Proportion of time at sea
  mutate(day_class = cut(perc_time_at_sea, seq(0, 100, 5),
    include.lowest = T
  )) %>%
  # by class of Proportion of time at sea
  group_by(day_class) %>%
  # calculate the proportion of daily benthic dives
  summarise(
    nb_benthic_class = sum(nb_daily_dives_type),
    nb_total_class = sum(nb_daily_dives),
    prop_class = nb_benthic_class / nb_total_class
  )

# plot
prop_at_sea <- prop_at_sea_dt %>%
  # intiate plot
  ggplot(
    .,
    aes(x = day_class, y = prop_class)
  ) +
  geom_bar(stat = "identity", fill = "#008c49") +
  guides(x = guide_axis(angle = 45)) +
  scale_y_continuous(
    limits = c(0, 1),
    # breaks=seq(0,0.6,0.1),
    labels = function(x) {
      scales::percent(abs(x), 1)
    }
  ) +
  annotate("text", x = 5.3, y = 0.4, label = "Departure") +
  annotate("text", x = 15, y = 0.5, label = "Arrival") +
  # geom_curve(
  #   data =   tibble(
  #     x1 = c(4, 17),
  #     x2 = c(1.5, 19.5),
  #     y1 = c(0.4, 0.5),
  #     y2 = c(0.31, 0.55)
  #   ),
  #   aes(
  #     x = x1,
  #     y = y1,
  #     xend = x2,
  #     yend = y2
  #   ),
  #   arrow = arrow(length = unit(0.08, "inch")),
  #   linewidth = 0.5,
  #   color = "gray20",
  #   curvature = 0.3
  # ) +
  geom_curve(
    data = tibble(
      x1 = c(4),
      x2 = c(1),
      y1 = c(0.4),
      y2 = c(0.33)
    ),
    aes(
      x = x1,
      y = y1,
      xend = x2,
      yend = y2
    ),
    arrow = arrow(length = unit(0.08, "inch")),
    linewidth = 0.5,
    color = "gray20",
    curvature = 0.3
  ) +
  geom_curve(
    data = tibble(
      x1 = c(16),
      x2 = c(19.7),
      y1 = c(0.5),
      y2 = c(0.57)
    ),
    aes(
      x = x1,
      y = y1,
      xend = x2,
      yend = y2
    ),
    arrow = arrow(length = unit(0.08, "inch")),
    linewidth = 0.5,
    color = "gray20",
    curvature = -0.56
  ) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.line = element_line(arrow = arrow(
      angle = 20,
      length = unit(.10, "inches"),
      type = "closed"
    ))
  ) +
  labs(
    x = "Time spent at sea (%)",
    y = "Proportion of benthic dives"
  )

Code

# display
prop_at_sea

Figure 5: Proportion of benthic dives performed during trip at sea expressed as a percentage of time spent at sea. Time spent at sea is classified into intervals of five percent for graphical representation.

Notes

• ?
• n = 422

Test to compare if the proportion of benthic dives at arrival is higher than the one at departure:

• \(H_{0}: p_{departure} \ge p_{arrival}\)
• \(H_{1}: p_{departure} \lt p_{arrival}\)

Code

# test to compare the proportion between arrival and departure
z_test <- prop.test(
  x = c(
    prop_at_sea_dt %>%
      filter(day_class == first(day_class)) %>%
      pull(nb_benthic_class),
    prop_at_sea_dt %>%
      filter(day_class == last(day_class)) %>%
      pull(nb_benthic_class)
  ),
  n = c(
    prop_at_sea_dt %>%
      filter(day_class == first(day_class)) %>%
      pull(nb_total_class),
    prop_at_sea_dt %>%
      filter(day_class == last(day_class)) %>%
      pull(nb_total_class)
  ),
  alternative = "less"
)
z_test

    2-sample test for equality of proportions with continuity correction

data:  c(prop_at_sea_dt %>% filter(day_class == first(day_class)) %>% pull(nb_benthic_class), prop_at_sea_dt %>% filter(day_class == last(day_class)) %>% pull(nb_benthic_class)) out of c(prop_at_sea_dt %>% filter(day_class == first(day_class)) %>% pull(nb_total_class), prop_at_sea_dt %>% filter(day_class == last(day_class)) %>% pull(nb_total_class))
X-squared = 5973, df = 1, p-value < 2.2e-16
alternative hypothesis: less
95 percent confidence interval:
 -1.000000 -0.236269
sample estimates:
   prop 1    prop 2 
0.3093339 0.5505404 

Code

# extract the z statistic
sqrt(z_test$statistic)

X-squared 
 77.28513 

\(H_{0}\) is rejected, so we can say that the porportion of benthic dives at arrival is significantly superior than departure!

Based on: http://www.sthda.com/english/wiki/comparing-proportions-in-r

and https://bookdown.org/danieljcarter/r4steph/z-test-for-proportions.html

4.7 Supplementary

4.7.1 Table 1

Code

# shark and orca area = NA if no location data
data_dive <- data_dive %>%
  mutate(
    shark_area = if_else(is.na(Lat), NA, shark_area),
    orca_area = if_else(is.na(Lat), NA, orca_area)
  ) %>%
  # add year of deployment
  group_by(id) %>%
  mutate(
    year = first(format(date, format = "%Y")),
    shark_and_orca = if_else(is.na(Lat),
      NA,
      if_else(shark_area + orca_area == 2, 1, 0)
    ),
    shark_or_orca = if_else(is.na(Lat),
      NA,
      if_else(shark_area + orca_area >= 1, 1, 0)
    )
  )

# build table
data_dive %>%
  group_by(year, id) %>%
  summarise(
    n_dives_ind_all = length(DiveNumber),
    n_dives_ind_w_loc = length(DiveNumber[!is.na(Lat)]),
    n_dives_ind_wo_loc = length(DiveNumber[is.na(Lat)]),
    n_benthic_dives_all = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic"])
    ),
    n_benthic_dives_w_loc = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" & !is.na(Lat)])
    ),
    n_benthic_dives_ind_in_shark = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        shark_area == 1 &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_in_orca = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        orca_area == 1 &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_shark_and_orca = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        shark_and_orca == 1 &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_shark_or_orca = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        shark_or_orca == 1 &
        !is.na(Lat)])
    ),
    hour_day_departure = ifelse(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      first(as.numeric(format(date_tz[DiveTypeName == "Benthic"],
        format = "%H"
      )))
    ),
    hour_day_arrival = ifelse(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      last(as.numeric(format(date_tz[DiveTypeName == "Benthic"],
        format = "%H"
      )))
    ),
    .groups = "drop"
  ) %>%
  group_by(year) %>%
  summarise(
    n_ind = length(id),
    n_dives_all = sum(n_dives_ind_all, na.rm = T),
    n_dives_w_loc = sum(n_dives_ind_w_loc, na.rm = T),
    n_dives_wo_loc = sum(n_dives_ind_wo_loc, na.rm = T),
    n_benthic_dives = sum(n_benthic_dives_all, na.rm = T),
    n_benthic_dives_proportion = sum(n_benthic_dives, na.rm = T) / sum(n_dives_all, na.rm = T),
    n_benthic_dives_in_shark_proportion = sum(n_benthic_dives_ind_in_shark, na.rm = T) / sum(n_benthic_dives_w_loc, na.rm = T),
    n_benthic_dives_in_orca_proportion = sum(n_benthic_dives_ind_in_orca, na.rm = T) / sum(n_benthic_dives_w_loc, na.rm = T),
    n_benthic_dives_in_shark_and_orca = sum(n_benthic_dives_ind_shark_and_orca, na.rm = T) / sum(n_benthic_dives_w_loc, na.rm = T),
    n_benthic_dives_in_shark_or_orca = sum(n_benthic_dives_ind_shark_or_orca, na.rm = T) / sum(n_benthic_dives_w_loc, na.rm = T),
    n_dives_mean = mean(n_dives_ind_all, na.rm = T),
    n_dives_plus_minus = "\u00b1",
    n_dives_sd = sd(n_dives_ind_all, na.rm = T),
    n_benthic_dives_mean = mean(n_benthic_dives_all, na.rm = T),
    n_benthic_dives_plus_minus = "\u00b1",
    n_benthic_dives_sd = sd(n_benthic_dives_all, na.rm = T),
    hour_day_departure_mean = psych::circadian.mean(hour_day_departure, na.rm = T),
    hour_day_departure_plus_minus = "\u00b1",
    hour_day_departure_sd = psych::circadian.sd(hour_day_departure, na.rm = T)$sd,
    hour_day_arrival_mean = psych::circadian.mean(hour_day_arrival, na.rm = T),
    hour_day_arrival_plus_minus = "\u00b1",
    hour_day_arrival_sd = psych::circadian.sd(hour_day_arrival, na.rm = T)$sd
  ) %>%
  gt() %>%
  tab_spanner(
    label = "All dives",
    columns = c(
      "n_dives_all",
      "n_dives_w_loc",
      "n_dives_wo_loc"
    )
  ) %>%
  tab_spanner(
    label = "Benthic dives",
    columns = c(
      "n_benthic_dives",
      "n_benthic_dives_proportion"
    )
  ) %>%
  tab_spanner(
    label = "% Benthic dives",
    columns = c(
      "n_benthic_dives_in_shark_proportion",
      "n_benthic_dives_in_orca_proportion",
      "n_benthic_dives_in_shark_or_orca",
      "n_benthic_dives_in_shark_and_orca"
    )
  ) %>%
  tab_spanner(
    label = "# dives",
    columns = c(
      "n_dives_mean",
      "n_dives_plus_minus",
      "n_dives_sd"
    )
  ) %>%
  tab_spanner(
    label = "# benthic dives",
    columns = c(
      "n_benthic_dives_mean",
      "n_benthic_dives_plus_minus",
      "n_benthic_dives_sd"
    )
  ) %>%
  tab_spanner(
    label = "Hour departure",
    columns = c(
      "hour_day_departure_mean",
      "hour_day_departure_plus_minus",
      "hour_day_departure_sd"
    )
  ) %>%
  tab_spanner(
    label = "Hour arrival",
    columns = c(
      "hour_day_arrival_mean",
      "hour_day_arrival_plus_minus",
      "hour_day_arrival_sd"
    )
  ) %>%
  tab_spanner(
    label = "per individual",
    columns = c(
      "n_dives_mean",
      "n_dives_plus_minus",
      "n_dives_sd",
      "n_benthic_dives_mean",
      "n_benthic_dives_plus_minus",
      "n_benthic_dives_sd",
      "hour_day_departure_mean",
      "hour_day_departure_plus_minus",
      "hour_day_departure_sd",
      "hour_day_arrival_mean",
      "hour_day_arrival_plus_minus",
      "hour_day_arrival_sd"
    )
  ) %>%
  # rename columns
  cols_label(
    year = "Year",
    n_ind = "# seals",
    n_dives_all = "Total",
    n_dives_w_loc = "w/ loc",
    n_dives_wo_loc = "w/o loc",
    n_benthic_dives = "Total",
    n_benthic_dives_proportion = "Proportion",
    n_benthic_dives_in_shark_proportion = "Shark area",
    n_benthic_dives_in_orca_proportion = "Orca area",
    n_benthic_dives_in_shark_or_orca = "Total",
    n_benthic_dives_in_shark_and_orca = "Overlap",
    n_dives_mean = md("Mean"),
    n_dives_plus_minus = "\u00b1",
    n_dives_sd = md("SD"),
    n_benthic_dives_mean = md("Mean"),
    n_benthic_dives_plus_minus = "\u00b1",
    n_benthic_dives_sd = md("SD"),
    hour_day_departure_mean = md("Mean"),
    hour_day_departure_plus_minus = "\u00b1",
    hour_day_departure_sd = md("SD"),
    hour_day_arrival_mean = md("Mean"),
    hour_day_arrival_plus_minus = "\u00b1",
    hour_day_arrival_sd = md("SD")
  ) %>%
  fmt_percent(
    columns = c(
      "n_benthic_dives_in_shark_proportion",
      "n_benthic_dives_in_orca_proportion",
      "n_benthic_dives_proportion",
      "n_benthic_dives_in_shark_or_orca",
      "n_benthic_dives_in_shark_and_orca"
    ),
    decimal = 1
  ) %>%
  fmt_number(
    columns = c(
      "n_dives_all",
      "n_dives_w_loc",
      "n_dives_wo_loc",
      "n_benthic_dives"
    ),
    decimal = 0
  ) %>%
  fmt_number(
    columns = c(
      "n_dives_mean",
      "n_dives_sd",
      "n_benthic_dives_mean",
      "n_benthic_dives_sd",
      "hour_day_departure_mean",
      "hour_day_departure_sd",
      "hour_day_arrival_mean",
      "hour_day_arrival_sd"
    ),
    decimal = 1,
    use_seps = FALSE
  ) %>%
  # alignement
  cols_align(
    columns = c(
      "n_dives_sd",
      "n_benthic_dives_sd",
      "hour_day_departure_sd",
      "hour_day_arrival_sd"
    ),
    align = "left"
  ) %>%
  # color rows
  opt_row_striping() %>%
  tab_options(table.width = pct(180))

											per individual
Year	# seals	All dives			Benthic dives		% Benthic dives				# dives			# benthic dives			Hour departure			Hour arrival
		Total	w/ loc	w/o loc	Total	Proportion	Shark area	Orca area	Total	Overlap	Mean	±	SD	Mean	±	SD	Mean	±	SD	Mean	±	SD
2004	26	280,155	215,452	64,703	9,291	3.3%	79.0%	100.0%	100.0%	79.0%	10775.2	±	3031.8	357.3	±	575.3	17.7	±	1.5	14.4	±	1.4
2005	35	281,636	226,566	55,070	16,413	5.8%	57.8%	100.0%	100.0%	57.8%	8046.7	±	3704.9	468.9	±	958.3	18.7	±	1.4	19.2	±	2.0
2006	29	219,773	181,667	38,106	5,219	2.4%	72.5%	100.0%	100.0%	72.5%	7578.4	±	3613.8	180.0	±	232.1	18.7	±	1.5	11.3	±	1.6
2007	30	235,691	235,691	0	5,829	2.5%	91.4%	100.0%	100.0%	91.4%	7856.4	±	3880.7	194.3	±	366.8	18.5	±	1.6	15.4	±	2.1
2008	35	252,268	239,800	12,468	3,460	1.4%	99.9%	100.0%	100.0%	99.9%	7207.7	±	3557.8	98.9	±	135.7	13.0	±	2.2	12.6	±	2.0
2009	24	149,429	135,010	14,419	3,646	2.4%	94.4%	100.0%	100.0%	94.4%	6226.2	±	3076.3	151.9	±	188.2	15.3	±	1.8	18.1	±	1.8
2010	35	248,956	248,956	0	3,702	1.5%	93.0%	100.0%	100.0%	93.0%	7113.0	±	3475.9	105.8	±	81.0	20.6	±	1.6	11.9	±	1.8
2011	32	244,228	232,301	11,927	3,883	1.6%	85.3%	100.0%	100.0%	85.3%	7632.1	±	3842.4	121.3	±	151.4	22.7	±	1.7	7.8	±	1.7
2012	33	256,295	256,295	0	4,233	1.7%	82.2%	100.0%	100.0%	82.2%	7766.5	±	4090.7	128.3	±	152.0	20.1	±	1.3	13.4	±	2.6
2013	31	265,526	265,526	0	5,149	1.9%	87.0%	100.0%	100.0%	87.0%	8565.4	±	4318.0	166.1	±	218.6	15.2	±	2.2	6.9	±	1.7
2014	22	141,911	141,911	0	6,679	4.7%	90.5%	100.0%	100.0%	90.5%	6450.5	±	3424.0	303.6	±	473.9	17.0	±	1.7	17.7	±	1.7
2015	30	226,090	226,090	0	4,380	1.9%	99.3%	100.0%	100.0%	99.3%	7536.3	±	3956.1	146.0	±	91.3	17.6	±	1.1	18.9	±	2.0
2016	26	209,870	209,870	0	4,003	1.9%	93.0%	100.0%	100.0%	93.0%	8071.9	±	4372.8	154.0	±	101.5	22.1	±	1.6	13.1	±	1.7
2017	10	75,326	75,326	0	1,648	2.2%	65.8%	100.0%	100.0%	65.8%	7532.6	±	3664.9	164.8	±	162.3	5.4	±	1.6	0.1	±	1.0
2018	17	146,244	141,809	4,435	2,733	1.9%	92.0%	100.0%	100.0%	92.0%	8602.6	±	3618.3	160.8	±	89.8	11.2	±	2.2	10.4	±	2.2
2019	7	32,823	32,823	0	769	2.3%	97.8%	100.0%	100.0%	97.8%	4689.0	±	1342.1	109.9	±	65.8	15.7	±	1.6	20.1	±	2.6

Notes

• lmk!

4.7.2 Figure 1

Let’s first download the map

Code

# with or without bathymetry
with_bathy <- F

# if bathy
if (with_bathy) {
  # check if background_ggoceanmaps exist
  if (file.exists("../export/background_ggoceanmap_zoom_out.rds")) {
    trip_zoom_out <- readRDS("../export/background_ggoceanmap_zoom_out.rds")
  } else {
    # using ggOceanMaps
    trip_zoom_out <- basemap(
      limits = c(170, -110, 30, 59),
      # limits = c(-155, -110, 30, 60),
      bathymetry = TRUE,
      shapefiles = "Arctic",
      rotate = TRUE,
      grid.col = NA
    )

    # Make the graticules:
    lims <- attributes(trip_zoom_out)$limits
    graticule <- sf::st_graticule(
      c(lims[1], lims[3], lims[2], lims[4]),
      crs = attributes(trip_zoom_out)$proj,
      lon = seq(-180, 180, 45),
      lat = seq(-90, 90, 10)
    )

    # Plot
    trip_zoom_out <- trip_zoom_out +
      geom_sf(data = graticule, color = "grey50")

    # reorder layers
    trip_zoom_out$layers <- trip_zoom_out$layers[c(1, 3, 2)]

    # save result
    saveRDS(trip_zoom_out, "../export/background_ggoceanmap_zoom_out.rds")
  }
} else {
  # using ggOceanMaps
  trip_zoom_out <- basemap(
    limits = c(170, -110, 20, 59),
    # limits = c(-155, -110, 30, 60),
    # land.col = 'black',
    # land.col = '#cdaa6d',
    bathymetry = FALSE,
    shapefiles = "Arctic",
    rotate = TRUE,
    grid.col = NA
  )

  # Make the graticules:
  lims <- attributes(trip_zoom_out)$limits
  graticule <- sf::st_graticule(
    c(lims[1], lims[3], lims[2], lims[4]),
    crs = attributes(trip_zoom_out)$proj,
    lon = seq(-180, 180, 45),
    lat = seq(-90, 90, 10)
  )

  # Plot
  trip_zoom_out <- trip_zoom_out +
    geom_sf(data = graticule, color = "grey70", linewidth = 0.3) +
    theme(
      panel.background = element_rect(fill = "#4292c6"),
      panel.ontop = FALSE
    )

  # reorder layers
  trip_zoom_out$layers <- trip_zoom_out$layers[c(2, 1)]
}

Tip

This part is mostly based on https://cran.r-project.org/web/packages/eks/vignettes/tidysf_kde.html. The author of the package confirms that his kernel density estimation is calculated using “meaningful” units such as UTM ¹, since the input is in a simple feature format.

Code

# data use to compute kernel density estimation
df_kernel_dens <- data_dive %>%
  # only with location data
  filter(!is.na(Lat)) %>%
  # select only the required columns
  select(lat = Lat, long = Long, id, DiveTypeName) %>%
  # create col to nicely display dive type
  mutate(origin = paste(DiveTypeName, "dives"))

Code

# transform data into sf object
df_kernel_dens_sf <- st_as_sf(
  df_kernel_dens,
  coords = c("long", "lat"),
  crs = st_crs(4326)
)

# make it's group_by origin
df_kernel_dens_sf <- group_by(df_kernel_dens_sf, origin)

# kernel density estimation
df_kernel_dens_sf_kde <- st_kde(df_kernel_dens_sf,
  H = diag(c(
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 1]
    ),
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 2]
    )
  ) / 4)^2
)

# https://github.com/r-spatial/sf/issues/1762
sf::sf_use_s2(FALSE)

# plot
to_plot <- trip_zoom_out +
  # rename axis
  labs(x = "Longitude", y = "Latitude") +
  # new scale
  new_scale_fill() +
  # kernel
  geom_sf(
    data = st_get_contour(
      # geospatial kernel
      df_kernel_dens_sf_kde,
      # probabilities
      cont = c(50, 80, 95, 99)
    ),
    # display
    aes(fill = label_percent(contlabel)),
    alpha = 0.7
  ) +
  # same colour bar
  scale_fill_viridis_d(option = "plasma") +
  # legend
  labs(fill = "Probs") +
  # no display alpha
  guides(alpha = "none") +
  # facet by origin
  facet_wrap(. ~ origin)

# reorder layers
to_plot$layers <- if (with_bathy) {
  to_plot$layers[c(1, 2, 4, 3)]
} else {
  to_plot$layers[c(1, 3, 2)]
}

# display
to_plot

Figure 6: Map of kernel density estimation representing the spatial distribution of Benthic, Drift, Forage, and Transit dives performed by northern elephant seals throughout their trip to sea. Probs represents different levels of probability density.

Notes

• lmk!

4.7.3 Figure 2

Code

# the last plot
prop_at_sea_trip <- data_dive %>%
  group_by(id, trip) %>%
  # time difference between the first date and the others (telling R to take diff b/w first date and all other dates)
  mutate(nb_days_departure = trunc(as.numeric(difftime(
    date, first(date),
    units = "days"
  )))) %>%
  # group by day/dive_type/seal
  group_by(nb_days_departure, DiveTypeName, id, trip) %>%
  # count number of dives
  summarise(nb_daily_dives_type = n(), .groups = "drop") %>%
  # group by day/seal
  group_by(nb_days_departure, id) %>%
  # count the total number of dives per day/seal
  mutate(nb_daily_dives = sum(nb_daily_dives_type)) %>%
  # order by seals/date
  arrange(id, nb_days_departure, trip) %>%
  # perform our calculation per seal
  group_by(id, trip) %>%
  # calculate of the Proportion of time since departure
  mutate(perc_time_at_sea = round(
    nb_days_departure / max(nb_days_departure) * 100, 1
  )) %>%
  # filter on benthic dives
  filter(DiveTypeName == "Benthic") %>%
  # add class of Proportion of time at sea
  mutate(day_class = cut(perc_time_at_sea, seq(0, 100, 5),
    include.lowest = T
  )) %>%
  # by class of Proportion of time at sea
  group_by(day_class, trip) %>%
  # calculate the proportion of daily benthic dives
  summarise(
    nb_benthic_class = sum(nb_daily_dives_type),
    nb_total_class = sum(nb_daily_dives),
    prop_class = nb_benthic_class / nb_total_class
  ) %>%
  # plot
  ggplot(
    .,
    aes(x = day_class, y = prop_class)
  ) +
  geom_bar(stat = "identity", fill = "#008c49") +
  guides(x = guide_axis(angle = 45)) +
  scale_y_continuous(
    limits = c(0, 1),
    # breaks=seq(0,0.6,0.1),
    labels = function(x) {
      scales::percent(abs(x), 1)
    }
  ) +
  facet_wrap(. ~ trip) +
  # annotate("text", x = 5.3, y = 0.4, label = "Departure") +
  # annotate("text", x = 15, y = 0.5, label = "Arrival") +
  # geom_curve(
  #   data =   tibble(
  #     x1 = c(4),
  #     x2 = c(1),
  #     y1 = c(0.4),
  #     y2 = c(0.33)
  #   ),
  #   aes(
  #     x = x1,
  #     y = y1,
  #     xend = x2,
  #     yend = y2
  #   ),
  #   arrow = arrow(length = unit(0.08, "inch")),
  #   linewidth = 0.5,
  #   color = "gray20",
  #   curvature = 0.3
  # ) +
  # geom_curve(
  #   data =   tibble(
  #     x1 = c(16),
  #     x2 = c(19.7),
  #     y1 = c(0.5),
  #     y2 = c(0.57)
  #   ),
  #   aes(
  #     x = x1,
  #     y = y1,
  #     xend = x2,
  #     yend = y2
  #   ),
  #   arrow = arrow(length = unit(0.08, "inch")),
  #   linewidth = 0.5,
  #   color = "gray20",
  #   curvature = -0.56
  # ) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.line = element_line(arrow = arrow(
      angle = 20,
      length = unit(.10, "inches"),
      type = "closed"
    ))
  ) +
  labs(
    x = "Time spent at sea (%)",
    y = "Proportion of benthic dives"
  )

Code

# display
prop_at_sea_trip

Figure 7: Proportion of benthic dives performed during trip at sea expressed as a percentage of time spent at sea during post breeding and post molting trip. Time spent at sea is classified into intervals of five percent for graphical representation.

4.8 Extra

Code

# import the csv export
data_plot <- readr::read_csv("../data/benthic_dive_extract.csv")

# plot
data_plot %>%
  ggplot(aes(x = time, y = CorrectedDepth)) +
  geom_path() +
  scale_y_reverse() +
  theme_minimal() +
  theme(
    axis.line.y = element_line(arrow = arrow(
      type = "closed",
      ends = "first",
      length = unit(0.1, "inches")
    )),
    axis.line.x = element_line(arrow = arrow(
      type = "closed",
      ends = "last",
      length = unit(0.1, "inches")
    ))
  ) +
  labs(
    x = "Time",
    y = "Depth (m)"
  ) +
  # geom_text(data = data.frame(
  #          xpos = median(data_plot$time),
  #          ypos = 30,
  #          hjustvar = 0,
  #          vjustvar = 0,
  #          annotateText = paste0("Corner index: ", round(unique(data_plot$corner_index),2),
  #                       "\n",
  #                       "Benthix index: ", round(unique(data_plot$benthic_index)),2)),
  #           aes(x=xpos,y=ypos,hjust=hjustvar,
  #               vjust=vjustvar,label=annotateText)) +
  scale_x_datetime(position = "top") +
  labs(subtitle = "ID: 2011034 - Dive: 29")

Footnotes

1. • https://gis.stackexchange.com/questions/64638/how-to-fill-in-the-parameters-for-kernel-density-estimation
   • https://digitalcommons.usf.edu/cgi/viewcontent.cgi?article=1020&context=siv
   • https://gis.stackexchange.com/questions/50692/choosing-projection-for-rasterization-of-random-latitude-and-longitude-data-in-n/50791#50791